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Abstract 

State-space models (SSMs) are increasingly used in ecology to model time-series 
such as animal movement paths and population dynamics. This type of hierarchical 
model is often structured to account for two levels of variability: biological stochastic- 
ity and measurement error. SSMs are flexible. They can model linear and nonlinear 
processes using a variety of statistical distributions. Recent ecological SSMs are of¬ 
ten complex, with a large number of parameters to estimate. Through a simulation 
study, we show that even simple linear Gaussian SSMs can suffer from parameter- 
and state-estimation problems. We demonstrate that these problems occur primarily 
when measurement error is larger than biological stochasticity, the condition that often 
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drives ecologists to use SSMs. Using an animal movement example, we show how these 
estimation problems can affect ecological inference. Biased parameter estimates of a 
SSM describing the movement of polar bears {Ursus maritimus) result in overestimat¬ 
ing their energy expenditure. We suggest potential solutions, but show that it often 
remains difficult to estimate parameters. While SSMs are powerful tools, they can 
give misleading results and we urge ecologists to assess whether the parameters can be 
estimated accurately before drawing ecological conclusions from their results. 


Keywords: Ocean Tracking Network, State-space model. Parameter estimability. Dynamic 
linear model 

1 Introduction 

State-space models (SSMs) are increasingly nsed in ecology and are becoming the favonred 


statistical framework for modelling animal movement and population dynamics (Buckland 


et ah, 2004 Patterson et akf 2008, McClintock et ah, 2012, Newman et al., 2014). SSMs 
are desirable because they are structured so as to differentiate between two distinct sources 
of variability: the biological or process variation (e.g., demographic stochasticity) and the 


measurement error associated with the sampling method (Patterson et ah, 2008 Newman 


et ah, 2014). Because marine observations are often associated with large measurement 


errors that can mask biological signals, much of the early development of SSMs in ecology 


was by marine ecologists and fisheries scientists (e.g., Newman, 1998 Jonsen et ah, 2003 


Sibert et ah, 2003). The SSM framework has since become a general approach to account 


for multiple levels of stochasticity when modelling time-series, making them increasingly 


popular in the terrestrial literature (e.g., Csillery et al., 2013 Fukasawa et ah, 2013 Flesch 


2014). Here, we demonstrate that even simple SSMs can be problematic. The model we 
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chose is often used to explain how SSMs can account for two levels of stochasticity (e.g., 


Newman et ah, 2014), yet, we show that it suffers from parameter- and state-estimation 


problems. 

SSMs are a type of hierarchical model, in which one level treats the underlying unobserved 
states as an autocorrelated process, while another level accounts for measurement error 


(Cressie et ah, 2009). The SSM framework is flexible, especially when htted with Monte 


Carlo methods such as particle hlters or Markov Chain Monte Carlo (MCMC). SSMs can 
be used to model a variety of linear and nonlinear processes, and can represent stochasticity 


with diverse statistical distributions (e.g., Pedersen et ah, 2011 McClintock et ah, 2012 


Albertsen et ah, 2015). The flexibility of the SSM approach allows ecologists to build complex 


models that describe the biological and measurement processes with levels of detail that were 
previously unattainable. 

While the SSM framework is flexible, much of its theoretical foundation is based on simple 


linear Gaussian SSMs (sometimes referred as normal dynamic linear models, see Newman 


et ah, 2014). An example of a simple univariate linear Gaussian SSM is the one we will use 


to demonstrate parameter-estimability problems: 


Measurement eq yt = Xt + €t , et rs-/ N{0, al) , where t > 1, > 0 (1) 

Process eq Xt = pxt-i -\- Pti Vt ^ N{0, cr^) , where t > 1, cr^ > 0 , —1 < p < 1, (2) 

where y = {yi,y2, ■■■,yn) are observed at regular time intervals t = ( 1 ,...,? 7 ,) for a 

time-series of length n and x = (xq, Xi,..., Xt,..., x„) are the true unobserved states, with 

Xq representing the initial state. An ecological example of such a time-series would be a 
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series of yearly population size estimates. For instance, Newman et al. (2014) use this 


model to introduce SSM for population dynamics with xt representing the true but unknown 
abundance of an animal population at time t, yt an unbiased observation of the population 
size at time t, and p the population growth rate. 

The origin of SSMs is intimately linked with the Kalman hlter, a recursive procedure to 
estimate the unobserved states based on inaccurate observations (e.g., estimating the true 
hsh abundance based on catch data). The Kalman hlter was developed to estimate states 


based on a model without unknown parameter values (Kalman, 1960). However, in eco¬ 


logical applications, most parameters need to be estimated (e.g., McClintock et ah, 2012). 


Fitting methods for SSMs, such as the Kalman hlter, are now used to facilitate both state 


and parameter estimation (Johnson et al., 2008). In many cases, SSMs are used to estimate 


variance parameters because they are designed to diherentiate measurement error from pro¬ 


cess stochasticity (Dennis et ah, 2006 Simmons et al., 2015). While estimating parameters 


is often a means to estimate the unobserved states (e.g., Johnson et al., 2008 Albertsen 


et ah, 2015), parameters themselves can be of interest because they describe the underly¬ 


ing dynamics of the system, or behaviour of the animal (e.g. Mills Flemming et al., 2010 


McClintock et ah, 2012). 


Estimability problems associated with SSMs and other hierarchical models have been dis¬ 


cussed in the population dynamics literature (e.g., Dennis et al., 2006 Knape, 2008). In 


particular, previous studies have emphasized how difficult it is to use SSMs to estimate den¬ 
sity dependence parameters ([Knape, 2008; Polansky et al., 2009) and to diherentiate process 


stochasticity from measurement error (e.g., Dennis et ah, 2006). However, the existence of 
parameter estimation problems have been largely overlooked in the movement literature, and 
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by those that use complex Bayesian SSMs. As SMMs are becoming the favoured framework 
for many ecological analyses (Buckland et ah, 2004 Patterson et ah, 2008t McClintock et al. 


2012 Newman et ah, 2014), and are gaining popularity in other helds (e.g., Cao et ah, 2014), 


it is timely to warn researchers of their weaknesses. 


Here, we use simulations to show that simple SSMs can have severe parameter-estimability 
problems that in turn affect state estimates. These problems are more frequent when the 
measurement error is large, the very condition under which SSMs are needed, and can 
persist even when we incorporate measurement error information. While our main estimation 
approach consists of maximizing the likelihood numerically through Template Model Builder 
(TMB; developed by Kasper Kristensen and available at www.tmb-project.org), we show that 
these problems persist across a wide range of platforms and statistical frameworks, including 
when the parameters and states are estimated via Bayesian methods. We use the polar bear 
{Ursus maritimus) movement data that led us to notice these problems to demonstrate the 
effect of estimation problems on the biological interpretation of results. Finally, we discuss 
techniques to diagnose and, when possible, alleviate estimability problems. 


2 Methods 

2.1 Demonstration of the problem 

When we £t models to data, we want the parameters to be identihable, which means that, 
given perfect data (e.g., an inhnitely long time-series), it is possible to learn the true values of 
parameters. Assessing parameter identihability is often difficult and a more attainable goal 
is to assess estimability. Estimability means that, given the data at hand, the method used 
to approximate the parameter yields a unique estimate. When the maximum value of the 
likelihood function occurs at more than one parameter value, the parameter is nonestimable. 
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The quality of parameter estimates can be assessed in terms of: its variance, measured 
over multiple repeated estimations; bias, the expected difference between the estimate and 
true value of the parameter; or mean square error, a composite of bias and variance. To 
demonstrate that the estimates of the parameters and states of SSMs can be inaccurate, we 
simulated a set of time-series using the model presented in eq. [T]j^ In all simulations, the 
values for the initial state, Xq, the measurement error, cie, and the correlation, p, were set to 
0, 0.1, and 0.7, respectively. In Appendix [^(Supplementary information), we explored other 
p values, including a simpler model where p is hxed to 1. Note that while this simpler model 


has fewer parameters to estimate, it is no longer stationary (Durbin and Koopman, 2001). 


To investigate whether the ratio of measurement to process stochasticity affected estimation, 
we simulated a range of cr^ values: (0.01, 0.02, 0.05,0.1, 0.2, 0.5,1). For each parameter set, 
we simulated 200 time-series each with 100 observations {n = 100). Analyses using longer 
time-series (n=500) are presented in Appendix [B] (Supplementary information). 

For each simulation, we estimated the parameters, 6 = (cXe, p, a^), and states, x, using the R 


R Core Team 

2015 

) package TMB. This R package is similar to AD Model Builder ( 

Fournier 


et al., 2012) in that it uses automatic differentiation and the Laplace approximation. Finding 


the Maximum Likelihood Estimate (MLE) of the parameters of a SSM requires the maxi¬ 


mization of the marginal distribution of the observations (Newman et ah, 2014). For the 
model presented in eq. [T]M this involves maximizing the following likelihood: 


Ma..p.a„|y) = /n p{yt\xt)p{,xt\xt-i)d^, 

^ t=l 
« n 

= / YlN{xt, a'^,)N{pxt-u al)d^. 


( 3 ) 

( 4 ) 


t=i 
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To get the marginal distribution, we integrate over the states, x = (xi,..., x„)'. In TMB, 
this integration is achieved using the Laplace approximation, which in turn also returns 


state estimates (Albertsen et ah, 2015). While we refer to state “estimation”, this pro¬ 


cess is sometimes called “prediction” because states can be interpreted as random variables 


(Newman et ah, 2014). In this example, we assumed that the initial state is known (i.e., 
Xq = 0), which should help the estimation process. In instances where the initial state value 


is unavailable, the initial state can be modelled as Xq ~ iV(/i, (Tq) (Durbin and Koopman 


2001). TMB calculates standard errors for the estimated parameters by using the inverse of 


the observed Fisher information, i.e. the Hessian of the log likelihood (similar to ADMB, 


see 


Fournier et ah, 2012). To calculate the 95% conhdence intervals (Cl), we multiplied the 
aforementioned standard errors by the 2.5 and 97.5*^ percentiles of the normal distribution 


(i.e., the quadratic approximation in Bolker, 2008). 


To demonstrate that the problem is widespread across different statistical platforms, we 
also htted the simulated data using two popular R packages: dim ([Petris 2010) and rjags 


(Plummer, 2014). dim uses the Kalman hlter for the state estimates and calculates the 


MLE with numerical optimization methods, rjags is an R interface to JAGS, a program 
that can be used to ht Bayesian hierarchical models using MCMC methods (Supplementary 
information: Appendix . 

We evaluated the parameter-estimation performance of SSMs by comparing the estimated 


and simulated values. Similar to Pedersen et ah (2011), we evaluated the state-estimation 


performance with the root mean square error (RMSE): 


RMSE = 


\ 


n 




( 5 ) 


t=i 
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where Xt is the estimated state at time t and Xt is the simulated (i.e., true) state at time t. To 
assess whether the state-estimation performance was affected by the parameter-estimation 
problems, we compare RMSE^, for which the parameters, 6 = p, a^), were also estimated, 
to RMSEg, for which the parameter values were hxed at the values used to simulate the 
data. 


To investigate the potential causes of the parameter-estimation problem, we explored the 
likelihood prohle for a subset of the problematic simulations. We used the same simula¬ 
tions and parameter values as above, with the exception that we only examined the most 
problematic values: = (0.01,0.02,0.05) (see Results). Because they are associated with 

high measurement error to process stochasticity ratios, these values also represent the con¬ 
ditions when SSMs are most needed. For each scenario (i.e., different values of cTj,), we 
randomly chose one simulation for which the RMSE^ was 50% larger than RMSE^. Again, 
we used TMB to estimate parameter values, 6 , and the states, x. To examine whether the 
estimation problems were associated with the simultaneous estimation of states and param¬ 
eters, we estimated parameters when the state values were hxed to their simulated values 
(Supplementary information: Appendix]^. As a hnal investigation of the causes of the es¬ 
timation problems, we show how these problems are associated with known limitations of 
the autoregressive-moving-average (ARMA) models (Supplementary information: Appendix 

0 . 

2.2 Incorporating measurement error information 


Many ecologists incorporate information on measurement error in their model by either 


hxing parameter values or, in a Bayesian framework, using informative priors (e.g. Jonsen 


et ah, 2003, 2005 Johnson et ah, 2008). We investigated whether hxing the measurement 












error resolved the parameter estimation problem. To do so, we fitted our simple likelihood 
(eqn. to the same simulations, but we hxed the standard deviation of the measurement 
equation to the value used to simulate the data, (7^ = 0.1. We only estimated the remaining 
parameters, 6^ = (p, o'??)- As above, we investigated the parameter estimates, RMSE of the 
states, and likelihood prohles. 

2.3 Ecological example 


The movement of many animals, such as birds, hsh and marine mammals, is a combination 
of the voluntary movement of the animal (active movement) and drift (passive displacement 
resulting from ocean or wind currents). Currents do not always direct animals towards 
their goals, and moving against currents may require a substantial amount of energy (e.g.. 


Weimerskirch et al., 2000). To understand how currents affect the behavioural strategies 


of an animal, it is necessary to distinguish between the voluntary movement of the animal 


and drift (Caspar et ah, 2006). The voluntary movement can then be used as a proxy of 


energy expenditure, or can be integrated into an energy budget model to assess the effects 
of movement on survival and reproduction (Caspar et ah, 2006 Molnar et al.| 2010, 2014). 
While developments in satellite telemetry are providing increasingly precise measurements of 
animal movement paths, it is difficult to differentiate between drift and voluntary movement 
because wind, ocean, and sea ice drift data are often associated with large errors (e.g.. 


Schwegmann et al. 

2011 

Fossette et al. 

2012) 


We noticed the estimation problems of linear Caussian SSMs when developing a model 
that would differentiate between the voluntary movement of polar bears and sea ice drift. 


Polar bears often move in the reverse direction of the sea ice drift (Mauritzen et al.l 2003 


Auger-Methe et al. 

2015 

) and sea ice drift can be associated with large errors ( 

Schwegmann 
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et al., 2011). As a proxy of energy expended by bears, we wanted to estimate the volnntary 


movement. As a first test, we developed a 2 dimensional SSM that accounts for error in ice 
drift data: 


Initial state 

xo ~ A^(^,Po) 


(6) 

Measurement eq 

yt = 

- Xt + St + et. 

et~iV(0,H) 

(7) 

Process eq 

Xt = 

= Txt_i + r/t. 

Vt ~ ^(o,Q), 

(8) 


where yt = is the measured daily displacement of the polar bear based on the GPS 

collar data, ] is the voluntary displacement of the polar bear, and ] is the 

daily sea ice drift experienced by the bear. Here, the measurement error, e^, is associated 
with the ice data, not the polar bear location data. The location data were determined by 


GPS, for which the error is negligible (< 30m, see Tomkiewicz et ah, 2010). For simplicity, 
we assumed that the two geographic coordinates are independent, thus: 


(9) 


Because eq. |6]j^ model displacements, the elements of H represent the measurement error in 
the sea ice drift data and those of Q are associated with the speed of the bear. Similar to 


Po = 

1 

o to 

0 

, H = 


0 

, Q = 

2 

^r),u 

0 

, T = 

1 

0 


0 

^0. 


0 



0 

2 


0 

pv 


7 in Jonsen et al. (2005), pu and represent the degree of autocorrelation in the random 


walk. To initialize the model we used p = [q] and ctq = 15. We chose 15 km as it is 
the standard deviation of the observed daily displacements of the polar bears in the u- and 
n-direction. 
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We used the daily movement of 15 polar bears collared in the Beaufort Sea in the spring 


of 2007-2011. The bears were immobilized with standard methods (Stirling et al., 1989) 
and equipped with Telonics Inc. (Mesa, AZ) collars. All capture and handling procedures 
were carried out in accordance with the protocols approved by the University of Alberta 
Animal Care and Use Committee for Biosciences. We used the Polar Pathhnder Daily 


25km Ease-Grid Sea Ice Motion Vectors (Fowler, 2003), which are daily estimates of sea ice 
displacements in the u- and u-directions of the Northern Hemisphere azimuthal equal-area 


EASE-Grid projection developed for polar sea ice data (Brodzik and Knowles, 2002). We 


used the same movement data and data handling procedures as in Auger-Methe et al. (2015), 
including interpolating the ice drift data at each bear location, assigning a drift value of zero 
for landfast ice, and excluding the three days after collaring to remove movements affected 
by handling. The only differences in the data used here, are that we excluded all bears that 
spent time on land and considered days with missing sea ice data as missing observations 
(i.e., we considered both y* and St as missing that day). 


Our goal was to use the SSM to estimate the energy expenditure of each bear. Our proxy 
was the total voluntary bear displacement: 


d = ^ \lxlu + xl^, ( 10 ) 

t=i 

where Xt,u and Xt^v are the estimates of the daily voluntary bear displacements in the u- 
and u-directions. The number of days, n, included in the time-series will affect our estimate 
of d. For consistency, we set n to be 342, the length of the shortest time-series across the 
15 bears. To assess the effects of estimation problems on our ecological interpretation, we 
simulated movement paths similar to those described by the polar bear data (Supplementary 
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information: Appendix . 

The code is available at https: //gitlab. oceantrack. org/otn-statistical-modelling-g;roup/ 
SSMestProblems and as Supplementary data. 

3 Results 

3.1 Simulations results 


According to the simulation results, parameter estimation was often inaccurate, and these 
problems affected the state estimates (Fig. [^. The parameter estimates were often far from 
their true values, and their distributions often bimodal (Fig. Supplementary Fig. A.l). 
In many cases, the estimates for cTe and p had peaks close to 0. The RMSE^ of the state 
estimates had either a bimodal distribution, or a long tail compared to that of the RMSE^ 


(Supplementary Fig. A.l). In other words, when the parameters were estimated, many 
replicates had much higher state estimate error than when the true parameter values were 
used (Fig. [^. In fact, 29.6% of the simulations had a RMSE^ value that was 50% larger than 
their RMSEg. When the simulations had high measurement error to process stochasticity 
ratios, the estimation problems for the states and two biologically relevant parameters, (p, cr^) 
were much higher (Fig. [^. The RMSE^ in some of these cases was close to 10 times greater 
than the simulated process stochasticity. 


Our supplementary analyses demonstrated that similar estimation problems occurred when 
dim and rjags were used (Supplementary information: Appendix]^. However, while the 
parameters estimated with rjags were often biased, their distributions did not contain a 
peak at 0. Increasing the length of the time-series improved parameter and state estimation 
(Supplementary information: Appendix 0. However, 500 time steps were insufficient to 
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completely eliminate problems. Our supplementary analyses also show that the problems 
are less apparent when p is close to 1, or when we used the simpler non-stationary local-level 
model, which hxes the value of p = 1 (Supplementary information: Appendix [A|). 

The likelihood prohles of a subset of the problematic simulations revealed that the likelihood 
was flat in some areas and sometimes bimodal or jagged (Fig. [^. The Cl of many param¬ 
eters excluded the true simulated value. Because the estimated measurement error of these 
simulations were close to 0, the estimated states were very close to the observations and far 
from their true simulated values (Fig. IP ,H,L). When the states were hxed to their simulated 
rather than estimated values, the likelihood prohles were unimodal and most Cl included 
the true parameter values, indicating that the problem lies in simultaneously estimating the 
states and the parameters (Supplementary information: Appendix [D|. 

3.2 Fixing the measurement error 

Fixing the standard deviation of the measurement error to the simulated value, = 0.1, 
helped reduce the estimation problems (Supplementary information: Appendix]^. RMSE^ 
values were much closer to RMSEg when the measurement error was hxed rather than esti¬ 
mated. In this case, only 5.0% of the simulations had a RMSE^ value that was 50% larger 
than their RMSEg. However, hxing the measurement error did not completely resolve the 
estimation problems. Some parameter estimates continued to be on the boundary of param¬ 
eter space and far from their simulated values. In addition, some likelihood prohles remained 
hat and some CIs spanned the entire parameter space (see Supplementary information: Ap¬ 
pendix]^ for more detail). 
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3.3 Ecological example 


The 15 polar bears studied used overlapping areas in the Beaufort Sea (Fig. [^), but 
their parameters estimates varied widely (Fig. In particular, three individuals had 

much lower estimated sea ice measurement error, with either < 0.01 and < 0.01. 
These three individuals had total voluntary displacement estimates that were on the higher 
end of the range (Fig. |B). These results are similar to those found when we simulated 
movement data similar to the real polar bear data (Supplementary information: Appendix 
1^. The simulations also showed that a few individuals would have < 0.01 and < 
0.01 and that these individuals would be associated with higher values of total voluntary 
displacement. 

4 Discussion 


Linear Gaussian SSMs, and approximations of them, are commonly used in the ecological 
literature to model animal movement (jJonsen et ah 2003 Johnson et ah, 2008; Patterson 


et ah, 2008) and population abundance (e.g., Wilson et ah, 2011 Flesch, 2014). These SSMs 
are often used to differentiate measurement error from process stochasticity and estimate the 


associated variance parameters (e.g., Sibert et ah, 2006 Wilson et ah, 2011 Flesch, 2014 


Albertsen et ah, 2015). Our results demonstrated that simple linear Gaussian SSMs can 


have severe parameter- and state-estimation problems, and that these problems can affect 
biological inferences. According to our simulations, estimation problems were more frequent 
when the measurement error was much larger than the process stochasticity. In such cases, 
the three estimated parameters were often far from their simulated values, which in turn 
resulted in inaccurate state estimates. The ARMA notation shows that when the measure¬ 
ment error is much greater than the process stochasticity there is parameter redundancy. 
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explaining why it is difficult to accurately estimate the parameters (Supplementary infor¬ 
mation: Appendix]^. Our simulations showed that hxing the measurement error to its true 
value helped, but did not completely solve the estimation problems, especially when the 
hxed measurement error was relatively large. This is particularly worrisome because SSMs 
are most needed when the measurement error is large compared to the process stochasticity, 
and this is the condition under which the largest estimation problems occur. 

The estimation problems are less critical when the measurement error is much smaller than 
the process stochasticity. While the measurement error estimates were often close to 0, the 
estimates for the other parameters, and those for the states, were generally accurate. As 
shown by the ARMA notation, when the measurement error is much smaller than the process 
stochasticity the model behaves as an AR(1) process, explaining why the measurement error 
estimates were often close to zero (Supplementary information: Appendix]^. In effect, the 
measurement error is ignored. However, when the measurement error is negligible compared 
to the process stochasticity, ignoring the effect of the measurement error is less likely to 
affect our interpretation of the biological process. 


Others have discussed estimation problems associated with htting simple linear Gaussian 
SSMs. A few recent ecological studies have reported difficulties when estimating variance 


parameters, including variance estimates close to 0 (Tittensor et ah, 2014 Simmons et al. 


2015). Dennis et al. (2006|, who transformed the stochastic Gompertz population model 


into a linear Gaussian SSM, noted that while the process stochasticity and measurement 
error parameters can be estimated, multimodal likelihood functions occur and can lead to 
erroneous estimates. They showed that the likelihood functions tended to have multiple 
peaks, including two peaks associated with either no process stochasticity or no measurement 
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error. While these two peaks can be local maxima, Dennis et ah (2006) noted that when 


there is substantial measurement error, one of these modes was often the global maximum. 


Knape (2008) extended the study of the Gompertz SSM to focus on the estimability of the 


density dependence parameter, an autocorrelation parameter similar to p. He found that 
the density dependence was generally not identihable in the presence of unknown process 
variability and measurement error, especially when the strength of the density dependence 
was close to 0. When the measurement error was known, the strength of density dependence 
was estimable but the estimates often remained biased. 

By extending the range of measurement error to process stochasticity ratios beyond those 


explored by Dennis et ah (2006) and Knape (2008), we demonstrate that relatively high 


measurement error can have dramatic effects on process parameter and state estimates, even 


when the measurement error is known. The results of Knape (2008) suggested that p val 


ues close to 0 would result in estimability problems (see also Forester et ah, 2007), which 


is not surprising. As the process becomes less autocorrelated it is harder to differentiate it 
from the temporally independent measurement error, suggesting that differentiating between 
measurement error and process stochasticity would require a large sample size when p is far 
from 1. However, our results demonstrated that estimation problems remained with rela¬ 
tively high autocorrelation, p = (0.7,0.99) and p hxed to 1, and relatively long time-series, 
n = (100,500) (see Supplementary information: Appendices [A|B| . These results emphasize 
that the parameters and states are only estimable for a narrow range of conditions. Both 
the analysis of the ARMA formulation of our SSM and our ecological example show that 
parameter estimability within linear Gaussian SSMs is a general issue, not one restricted to 
the stochastic Gompertz population model. In fact, these problems extend to some nonlinear 
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SSMs. For example, some of the estimated parameters of the nonlinear population SSMs of 


de Valpine and Hastings (2002) had considerable bias when measurement error was large rel¬ 


ative to process variability, de Valpine and Hilborn (2005) showed that their advance Monte 


Carlo kernel likelihood method could not differentiate between the process stochasticity and 


measurement error of the nonlinear Schaefer population model, and Polansky et ah (2009) 
found similar problems in the theta-Ricker model. 

Left undiagnosed, biased parameter estimates will mislead conclusions based on the prob¬ 
lematic model parameters and may affect our interpretation of the other model parameters. 


the state estimates, and other derived values (Cressie et ah, 2009 Lele, 2010). For example. 


stochastic population SSMs with negatively biased estimates of the process stochasticity will 


underestimate extinction risk (Lindley, 2003). In our polar bear example, erroneous esti¬ 


mates of measurement error and process stochasticity biased the state estimates and proxy 
for energy expenditure. Thus, even if the parameter values per se are not of interest, esti¬ 
mation problems need to be diagnosed because their effect on state estimates are likely to 
affect results of ecological importance. 

The hrst step to avoid these biased inferences is to detect the potential for parameter estima- 
bility problems, which can be done through a variety of practical means. Our simulations 
demonstrated that estimates at the boundary of parameter space can be indicative of a 
problem. For our polar bear example, we detected the estimation problem because we had 
no reason to believe that the three bears with sea ice measurement error close to 0 used 
different sea ice than the other bears. These three bears were exposed to similar levels of 
sea ice drift as other bears and were not geographically or temporally isolated from them. 


Investigating the likelihood prohle can also help detect estimation problems (Dennis et al. 
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2006, 2010 Ives et al., 2010). Indeed, the likelihood prohles of our problematic simulations 


had flat sections and multiple modes. However, in a Bayesian framework, the estimation 
problems can be obscured by the use of vague priors, as these can smooth the likelihood 


and affect inference (Bindley, 2003 Dennis et al., 2006; Lele and Dennis, 2009 Lele, 2010). 
When we used JAGS to estimates parameters, we had no estimates at the boundary and the 
posterior distributions of most parameters were unimodal, and yet, the estimates were biased 
(Supplementary information: Appendix]^. A useful way to evaluate the model’s capacity 
to separate process and measurement error parameters, is to assess the extent of correlation 
between these estimates (see Supplementary information; Appendix for details). In the 
maximum likelihood context, a plot of the likelihood surface can reveal a correlation pattern 


symptomatic of an identihability issue (de Valpine and Hilborn, 2005 Polansky et ah, 2009). 
In a Bayesian context, a plot of the joint posterior samples of these two parameters can reveal 
similar correlation patterns (Supplementary information: Appendix]^. While few methods 


have been developed to formally assess parameter identihability problems, data cloning (Lele 


et ah, 2010; Campbell and Lele, 2014) and the symbolic method (Cole, 2012 Newman et al. 


2014) are promising avenues. 


How can we avoid these estimability problems? In many cases, a larger sample size can 
help (see Supplementary information: Appendix |^. In particular, Dennis et al. (2010) 
demonstrated that sampling replicates can substantially improve the capacity of SSMs to 
differentiate process stochasticity from measurement error, and that it may be advantageous 
to design monitoring programs with multiple replicate counts per survey rather than increas¬ 
ing the length of the time series (i.e., number of times the survey is conducted). However, 
for many observational studies, ecologists are limited in their ability to gather more data 


Dennis et al. (2010 
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and, for movement data, it is often impossible to have replicates of location estimates. An 
alternative is to incorporate information on the measnrement error. As we demonstrated in 
onr simnlation stndy, when we £x the measnrement error to its trne valne, the estimates of 


the other parameters improved. While some parameter-estimation problems persisted, their 


effect on the state estimates diminished snbstantially. Similarly, de Valpine and Hilborn 


(2005) demonstrated that knowing the ratio of process to measurement variance would im¬ 
prove parameter estimates. In a Bayesian framework, specifying informative priors for the 
measurement error could help make the other parameters identihable and improve the state 


estimates (Bindley, 2003 Cressie et ah, 2009, but see Lele and Dennis 2009). Another al¬ 
ternative is to estimate the measurement error and process stochasticity outside of the SSM 
framework using the principle that the measurement error is uncorrelated over time whereas 


the process stochasticity is temporally correlated (Dowd and Joy, 2011). Estimating the 
measurement and process standard deviations offline reduces the number of parameters to 
estimate within the SSM framework. Using restricted maximum-likelihood, which treats 
hxed-effects parameters (e.g., p) and variance components (e.g., cx^) differently, can also 


be valuable to remove bias in SSM estimates (Dennis et ah, 2010). When the estimation prob¬ 
lem results in variance estimate close to 0, one can limit the estimate to interior (non-zero) 


solutions (Dennis et al., 2006 Knape, 2008). In particular, Dennis et ah (2006) suggested 
trying a variety of starting values for the optimizer used to numerically maximize the like¬ 
lihood and eliminating all solutions that involve variance with near 0 values, even if one of 
these is the global maximum. Finally, restructuring the model can help reduce the problem. 
For example, in the polar bear example, we could create a population model with a single 
measurement error parameter for all bears. Even if the process variability continues to differ 
between individuals, using one measurement error term for all bears signihcantly decreases 
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the number of parameters to estimate and increases the amount of data with which the mea¬ 
surement error term is estimated. As a general rule decreasing the number of parameters to 
estimate and increasing the amount of data will help reduce estimability problems. 


Not all parameters are equally affected by estimation problems. Forester et ah (2007), 
who developed a linear Gaussian SSM for animal movement, demonstrated that coefficient 
parameters associated with covariates and an intercept in the measurement equation are 
easier to separate than process autocorrelation (equivalent to p), measurement error and 
process stochasticity. Note, however, that all of these parameters had cases associated with 
estimation problems. For example, the coefficient estimates were biased when their true 


simulated value was not equal to zero. Humbert et ah (2009) suggested that in the case 
of exponential growth SSMs the population trend parameter, similar to an intercept in 
the process equation, was often well estimated and that increasing the precision of the 
abundance estimates and the length of the time series, more than the completeness of the 
time series, could increase the performance of the SSM. This further indicates that ecologists 
should closely consider model formulation, and that the estimability of parameter should be 
assessed. 


If we cannot resolve the parameter estimation problem, we need to account for its potential 
effect on our inference. One way to account for the estimation uncertainty is to use a 


parametric bootstrap to get CIs on the parameter and state estimates (Dennis et al., 2006 


Forester et al., 2007). These bootstrap CIs require simulating the model using the estimated 
parameter values and re-£tting the model to each simulation. The 2.5**^ and 97.5**^ quantiles 
of the estimated parameters and states then becomes the 95% Cl. These CIs differ from 
those we calculated from the standard deviation reported by TMB. However, because TMB 
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is orders of magnitude faster than MCMC methods (Albertsen et ah, 2015), implementing 


these parametric bootstrap CIs would be computationally feasible, even for complex models. 
Note, however, that the variability in the estimates of our simulations suggests that these 
CIs would be large and would often approach the boundary of parameter space. 

5 Conclusion 

We demonstrated that even simple linear Gaussian SSMs can have parameter estimability 
problems and that these problems can affect our ecological interpretation. As parameter 
estimability problems have been observed in other hierarchical models and because the ratio 
of information content to model complexity is expected to decrease with increasing numbers 


of hierarchies (Lele and Dennis, 2009 Lele, 2010), it is likely that these problems could occur 


in more complex forms of SSMs. Estimating individual variance components is notoriously 
difficult. SSMs do not escape this difficulty. While estimability problems have been discussed 


in the context of a few specihc population dynamics SSMs (e.g., Dennis et ah, 2006; Knape 


2008 Polansky et al., 2009), the voluminous literature on SSMs has paid relatively little 


attention to these problems. Such limited appreciation of the estimation problem is particu¬ 
larly dangerous because SSMs are usually advertised as providing the means to differentiate 


process from measurement variability (e.g., de Valpine and Hastings, 2002 Patterson et al. 


2008 Ahrestani et ah, 2013) 


It is timely to warn ecologists of these difficulties. SSMs are becoming the favoured frame¬ 
work for animal movement and population dynamics. SSMs used in ecology are becoming 


increasingly complex (e.g., McClintock et al., 2012). In addition, tools to apply SSMs to 
data are becoming increasingly available. For example, R now provides a variety of packages 


that £t SSMs (Petris and Petrone, 2011). Until recently, SSMs were applied by statisti- 
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cians or by ecologists with a strong statistical background. These researchers were more 
likely to be aware of potential estimability problems than most ecologists. Researchers have 
questioned whether ecologists have sufficient statistical training to properly implement hier¬ 
archical models and have suggested that universities should start including advanced courses 
in statistical modelling in their ecological programs (e.g., Dennis et ah, 2006 Cam, 2012). 
If the limitation of SSMs are not emphasized, the better accessibility of tools to £t these 
increasingly complex models are likely to lead to many undiagnosed estimation problems 
and incorrect conclusions. 

While SSMs are powerful tools, they can give misleading results if they are misused. We 
believe it is important for ecologists to be aware of the potential estimation problems of SSMs. 
Investigating the likelihood prohle, incorporating information on measurement error, and 
accounting for estimability uncertainty are all good hrst steps. However, we urge statisticians 
to develop further tools that can be used to diagnosed such problems and these should be 
readily available along with the tools to £t SSMs. 
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Figure Changes in parameter estimates and state RMSEs associated with varying the 
measnrement error to process stochasticity ratios (cTe/a^) in the simnlations. A-C) The 
boxplots represent the distribntion of the parameter estimates (a^, p, a^) and the pink circles 
represent the trne (simnlated) valnes. D) The grey boxplots represent the distribntion of the 
RMSE of the model htted nsing the estimated parameter valnes, while the pink boxplots 
represent the RMSE when the model is htted nsing the trne parameter valnes. 

Figure Log likelihood prohles for problematic simulations. In the hrst three columns, the 
curve represents the log likelihood when the focal parameter is hxed (the other parameter 
are optimise to maximise the log likelihood). The dash lines are the true parameter values 
(i.e., value used for the simulation), the full lines are the maximum likelihood estimates and 
the grey bands represent the 95% CL The last column shows the time-series. The black 
lines represent the observations, yt, fhe red lines the simulated true states, Xt, and the grey 
dashed lines the estimated states, Xt. 


Figure Polar bear movement, parameter estimates of the polar bear sea ice model, and 
estimates of the total voluntary bear displacement. A) Locations of the 15 polar bears used 
in the analysis, with colours representing different individuals. The map was created in 


R (R Core Team, 2015) using the Northern Hemisphere azimuthal equal-area EASE-Grid 


projection developed for polar sea ice data Brodzik and Knowles (2002). B) Estimated total 


voluntary displacement over 342 days. C-H) Parameter estimates of the polar bear sea ice 
models. The different colours in panels B-C represent the three individuals for which either 


< 0.01 or m,, < 0.01. 
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A Some problems persist when p is close to 1 and in the local-level model 


Knape (2008), who investigated a linear Gaussian SSM describing the stochastic Gompertz 


population model, found that an autocorrelation parameter similar to p was harder to esti¬ 
mate when its true simulated value was close to 0. p appears to be problematic especially 


when IpI < 0.5. Similarly, Forester et ah (2007), who developed a linear Gaussian SSM for 
animal movement, noticed that when the autocorrelation was close to 1 (i.e., 0.95) there was 
less estimation problems than when it was close to 0 (i.e., 0 or 0.2). This is not surprising. 
As the process becomes less autocorrelated it is harder to differentiate it from the temporally 
independent measurement error. As such, we focussed on investigating whether the estima¬ 
tion problems remained when the autocorrelation parameter was relatively high. In the main 
text, we have presented the results when p = 0.7. In this Appendix, we investigated higher 


p values. First, we recreated the same simulation study as in section [2T] of the main text, 
except that we used p = 0.99 in the simulations. Second, we used a simpler model called 


the local model, which is sometimes referred as the random walk plus noise (e.g., Petris 


2010 ): 


Measurement eq yt = Xt + et, 


et ~ A^(0, al ), where t > 1, > 0 (11) 


Process eq Xt = Xt-i + Pt, Pt ~ a'i ), where t > 1, > 0. (12) 


The only difference with the model presented in eq. [I]j^ is that there is no p parameters, 
which is the equivalent of hxing p = 1. Note that while this simpler model has fewer 


parameters to estimate, it is no longer stationary (Durbin and Koopman, 2001). Following 


the methods described in section 2.1, we simulated and htted this simpler model. 
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The parameter and state estimates improved in simnlations with p = 0.99 (Fig. A.2). 
In particular, the estimates for the process parameters (p, a,,) were much closer to their 
simulated values than when the simulated p value was 0.7 (e.g., compare Fig. A.2p -C 
to Fig. A.lp -C). In some cases, this translated in better state estimates (e.g., compare 
Fig. A.2p to Fig. A.lj d). Similarly, using the local model improved parameter and states 


estimates (Fig. A.3). However, this did not completely eliminate the estimation problems. 
When the measurement error was much larger than the process stochasticity, = 10 
some of the parameter estimates remained on the boundary of parameter space (e.g.. Fig. 
A.2|4 and Fig. A.3p). In the case of simulations with p = 0.99, the estimated p was close to 


0 (Fig. A.2p,F) and some state estimates were far from those estimated when the parameter 


values were known (e.g. Fig. A.2D,H). 
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Figure A.l: Histograms of the parameter estimates and of the RMSE of the state 
estimates when p = 0.7. This is a more detailed visualization of the results 
presented in Fig. [T] of the main text. Each row represents the results of 200 
simulations for a set of parameter values. For the hrst three columns, the vertical 
lines represent the parameter values used in the simulations, with black lines 
used for values that remained constant, = 0.1 and p = 0.7, and red lines for 
values that changed between sets, cr^ = (0.01,0.02,0.05,0.1,0.2,0.5,1). In the 
last column, the grey histograms represent the RMSE of the model htted using 
the estimated parameter values, while the blue histograms represent the RMSE 
when the model is htted using the true values. 
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Figure A.2: Histograms of the parameter estimates and of the RMSE of the states 
when p = 0.99. Each row represents the results of 200 simulations for a set of 
parameter values. For the hrst three columns, the vertical lines represent the 
parameter values used in the simulations, with black lines used for the values 
that remain constant for all simulation sets, Ug = 0.1 and p = 0.99, and red lines 
for values that change between set, cx^ = (0.01,0.02,0.05,0.1,0.2,0.5,1). In the 
last column, the grey histograms represent the RMSE of the model fitted using 
the estimated parameter values, while the blue histograms represent the RMSE 
when the model is fitted using the simulation values. 
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Figure A.3: Histograms of the parameter estimates and of the RMSE of the states 
for a set of simulations of the local-level model. Each row represents the results 
of 200 simulations for a set of parameter values. For the hrst three columns, the 
vertical lines represent the parameter values used in the simulations, with black 
lines used for the values that remain constant for all simulation sets, = 0.1, and 
red lines for values that change between set, = (0.01, 0.02,0.05, 0.1, 0.2, 0.5,1). 
In the last column the grey histograms represent the RMSE of the model fitted 
using the estimated parameter values, while the blue histograms represent the 
RMSE when the model is fitted using the simulation values. 
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B Some problems persist with longer time-series 


Using longer time-series can considerably reduce estimability problems. In this appendix, 
we investigate whether the estimation problem persisted with longer time-series. We wanted 
to use a length of time-series that is relevant to real ecological examples, keeping in mind 
that, generally, population abundance time-series are shorter than movement time-series. 
For example, a recent population study using bird count data was limited to 45 to 52 counts 


(Simmons et ah, 2015). As such previous simulation studies for population dynamics SSMs 


limited their time-series to 30 and 100 time steps (Dennis et al., 2006 Knape, 2008 Humbert 


et ah, 2009). In fact, Dennis et al. (2006) mentioned that times-series of 100 time steps were 


unrealistic for ecological data of the type and thus our example in the main text (n = 100) 
is likely underestimating the frequency of estimation problems. Movement time-series are 


generally of longer length. For example, the simulation study of Forester et al. (2007) used 


350 steps and their real movement data ranged from 265-390 locations. Our polar bear 
time-series ranged from 342 to 365. In addition, with technological advancements movement 
time-series are becoming much longer. Thus, to look at time-series more representative of 


movement time-series we conducted the same simulation analysis as in section 2.1 of the main 
text, with the only difference being that our time-series have 500 observations (n = 500) 
rather than 100 {n = 100). 

When time-series of 500 steps were used, the estimation problems were reduced (compare Fig. 


B.l to Fig. A.l). In particular, when Ue < 2 cr^, we had fewer cXe estimates at the boundary 


of parameter space (e.g., compare Fig. to Fig. A. Ip ) and the estimates of and p 

are closer to the simulated values (e.g., compare Fig. B.1T,0 to Fig. A.1[ T,0). However, 
when the measurement error is large compared to the process stochasticity, > 5 (7,,, many 
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(Te estimate remained at the boundary of parameter space (Fig. B.14,E), many cr^ estimates 
were positively biased (Fig. B.lp,G), and many p estimates were negatively biased, with 


values close to 0 (Fig. B.1B,F). In addition, when the parameters were estimated, many 
replicates had higher state estimate error than when the true parameter values were used 
(Fig, |b 3) ,H), indicating that biases in parameter estimates continued to affect the state 
estimates. Overall, these results suggest that longer time-series do improve the estimability of 
some parameters and states, but that to have reliable estimates when the measurement error 
is much larger than the process stochasticity would require much longer time-series. 
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Figure B.l: Histograms of the parameter estimates and of the RMSE of the states 
for a set of simulations with 500 time steps. Each row represents the results of 
200 simulations for a set of parameter values. For the hrst three columns, the 
vertical lines represent the parameter values used in the simulations, with black 
lines used for the values that remain constant, = 0.1 and p = 0.7, and red lines 
for values that change between set, cx^ = (0.01,0.02,0.05,0.1,0.2,0.5,1). In the 
last column, the grey histograms represent the RMSE of the model htted using 
the estimated parameter values, while the blue histograms represent the RMSE 
when the model is htted using the simulation values. 
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C Results with other R packages 


We chose to use TMB for most analyses because it is a fast and flexible package that can be 


used to fit a variety SSMs to data (e.g., Albertsen et al-f 2015). To verify that the estimation 
problems are not limited to this package and are general problems associated with linear 
Gaussian SSMs, we also used two additional packages to reproduce the simulation study 


explained in section 2.1 of the main text. First, we used dim (Petris, 2010), a package that 


was used in recent ecological studies to £t SSMs to data (Tittensor et ah, 2014 Simmons 


et ah, 2015). The package dim maximises the log likelihood numerically and has functions 


to estimate the states via Kalman filter and smoother. To be consistent with TMB, we 


used the Kalman smoother, which takes into account all observations (see Albertsen et al. 


2015). 


Second, we used rjags (Plummer, 2014), which is an R interface to JAGS, a program that 


allows for the analysis of Bayesian hierarchical models using Markov Ghain Monte Garlo 
(MGMG) methods. Unlike TMB and dim, rjags requires the specification of priors for the es¬ 
timated parameter. We used the vague priors: HaifN(0, cr^ = 10000), p ~ Uniform(0,1), 

and cr^ ~ HalfN(0, = 10000). We used two chains, each with 50 000 adaptation steps, 50 
000 burn in steps and 50 000 saved steps. For each chain we kept 1 every 500 steps. 

The resuits when we used dim are neariy identicai to those when we used TMB (compare Fig. 


G.l to Fig. A.l). The oniy difference, is that when cxe > 5 cr,, a few iess repiicates had 


(Te values close to 0, p values close to 0, and positively biased ar, values. However, these 
differences were small and some were accompanied by other biases, such as more p values 
close to -1. Overall, the conclusion made for TMB in the main text hold true for dim. 
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In contrast, the results from rjags are different from TMB. In particular, the distribution 
of estimates were unimodal and there was no estimates close to the boundary of parameter 
space (i.e., no close to 0). However, the peak of estimates was often far from the simulated 


value (Fig. C.3), indicating that both parameter and state estimates were often biased. This 
is potentially due to the fact that vague priors can influence the results and smooth the 


peculiarities of the likelihood (Dennis et ah, 2006). This effect could explain why estimation 
problems have been less easily detected in recent SSMs studies, which often uses complex 
Bayesian SSMs. In addition, the posterior distributions were more unimodal and not as 


flat as the likelihood prohles produced by likelihood-based methods (compare Fig. C.3 to 
Fig. §. Note that these results exclude the 37 replicates out of 1400 simulations that did 
not converge (scale reduction factor of any parameter > 1.1), and thus would have been 
deemed problematic with such metrics. Overall, the results from JAGS indicate that using 
Bayesian methods does not £x the estimation problems of the linear Gaussian SSMs, and in 
fact might have made them harder to detect. See Appendix]^ for a more detailed discussion 
of diagnostic tools. 
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Figure C.l: Histograms of the parameter estimates and of the RMSE of the states 
when dim is used to £t our SSM to a set of simulations. Each row represents the 
results of 200 simulations for a set of parameter values. For the hrst three columns, 
the vertical lines represent the parameter values used in the simulations, with the 
black lines used for the values that remain constant, Ue = 0.1 and p = 0.7, and the 
red lines for values that change between set, = (0.01, 0.02,0.05, 0.1, 0.2, 0.5,1). 
In the last column, the grey histograms represent the RMSE of the model htted 
using the estimated parameter values, while the blue histograms represent the 
RMSE when the model is htted using the simulation values. Note that the pa¬ 
rameters of 2 out of the 1400 simulations could not be estimated due to singularity 
problems. 
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Figure C.2: Histograms of the parameter estimates and of the RMSE of the states 
when rjags is used to £t our SSM to a set of simulations. Each row represents 
the results of 200 simulations for a set of parameter values. For the hrst three 
columns, the vertical lines represent the parameter values used in the simulations. 
The black lines are used for the values that remain constant for all simulation 
sets: (Tg = 0.1 and p = 0.7. The red line for that values that change in between 
set: = (0.01,0.02, 0.05, 0.1, 0.2, 0.5,1). In the last column the grey histograms 

represent the RMSE of the model htted using the estimated parameter values, 
while the blue histograms represent the RMSE when the model is htted using the 
simulation values. Note that 37 out of 1400 simulations did not converged (i.e. 
the potential scale reduction factor was >1.1 for one of the parameters). 
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Density Density Density 



Figure C.3: Posterior distribution for problematic simulations. The first three 
columns, the curve represents the posterior distribution for the estimated pa¬ 
rameters. The dash lines are the true parameter values (i.e., value used for the 
simulation), the full lines are the mean value. The last column shows the time- 
series. The black lines represent the observations, i/t, the red lines the simulated 
true states, Xt, and the grey dashed lines the estimated states, x*. Note that 
these three simulations converged (the potential scale reduction factor was <1.1 
for all parameters). 
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D Estimating parameters when we know the true states 


In this Appendix we investigate whether the parameter-estimation problems are associated 
with estimating the parameters at the same time as the states. To do so, we estimated the 
parameters when the states were fixed to their true simulated values. We focused on the 
same problematic simulations as those explored in the main text. When the states were 
known, the parameter estimates were close to the simulated values and the CIs included the 


true simulated values (Fig. D.l). In addition, the likelihood profiles were unimodal, demon¬ 
strating the type of likelihood profiles one would expect for well-behaved models. These 
results suggest that the problems lie in estimating both the states and the parameters. 
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Figure D.l: Log likelihood profile for problematic simulations when the state 
values are hxed to the simulated values. In the hrst three columns, the curve rep¬ 
resents the log likelihood when the focal parameter is hxed (the other parameter 
are optimize to maximize the log likelihood). The dash lines are the true param¬ 
eter values (i.e., value used for the simulation), the full lines are the maximum 
likelihood estimates and the grey bands represent the 95% CIs. The last column 
shows the time-series. The black lines represent the observations, yt, and the red 
lines the simulated true states, Xt- Note that these are the same simulations as 
in Fig. 
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E Reformulating our simple SSMs with the ARMA notation 

For certain parameter values, the SSM can behave either as a white noise process or as an 
AR(1) process. In both cases, the SSM formulation of the model will be over-parameterized, 
and will lead to estimation problems. To see this, we can rewrite our SSM as an ARMA(1,1) 
model. First, we can combine eq. and and reparametrize the model in terms of ep 

Vt = pxt-i + rjt + et, et ~ ~ A(0,a^) (13) 

= p{xt-i + Vt-i - Vt-i) +Pt + €-t (14) 

= pyt-i -h pet-1 + Pt + et (15) 


Since pt and e* are independent and normally distributed, their sum, ut = pt + et, follows a 
normal distribution with zero mean and variance = cr^+cr^. Now, if we let Ut-i ~ N{0, cr^), 
we can rescale its variance such that: 


D 


yt = pyt-i + P 


cr. 




zUt-i Ut. 


(16) 


This is an ARMA(1,1) process with AR parameter 0 = p, MA parameter 'ip = p 
variance parameter cx^. 


\AI+F 


, and 


When (Te (Tj^, then ip is small. Thus the process behaves as an AR(1) process with 


parameters p and cr^. This is the case for our simulations with cr,, = 0.5,1 (Fig. A.l). When 


(Te ^ CT„, then 'ip K, (p and hence there is parameter redundancy in the model (Box et al. 


2008). In this case, the process closely resembles white noise. This is the case for our most 


problematic simulations (a,, = 0.01,0.02,0.05, see Figs. A.l 2). 
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F Polar bear and sea ice simulations 


To demonstrate that the polar bear and sea ice model has estimation problems and show 
how these problems may affect our interpretation of our proxy of energy expenditure, we 
simulated movement using model described in eq. [6]j^ We simulated 500 movement paths 
with n = 342 using the parameters estimated with the polar bear data (Table [FT| . We also 
used the initial state values estimated with the polar bear data: xq = [-I'.os]- We simulated 
the sea ice displacement in the u- and u-direction using normal distributions with mean and 
standard deviation values based on the sea ice drift experienced by the polar bears in our 
sample: St,u ~ = 1.49, a = 4.97) and St,u ~ 77(p = 2.26, a = 3.97). As for the empirical 


data, we estimated parameters and the total voluntary displacement (see eq. 10). 


Table F.l: Parameter estimates for the polar bear sea ice empirical data. These 
are the parameter values used in the simulations. 


Parameters 

mean 

sd 


5.53 

3.10 


5.79 

2.66 

Pu 

0.635 

0.128 

Pv 

0.685 

0.123 


9.66 

2.47 

^7J,V 

8.56 

2.33 


Our simulation results show that the model appears to have similar estimation problems 
than the simple SSM extensively studied in the manuscript. In particular, we have a few 


simulations where the estimates of and are close to 0 (Fig. F.l), something that 


was also noticeable in the empirical data (Fig. [^. Parameter estimates close to 0 can be 
associated with estimation problems, which in turn can affect the state estimate and thus 


the estimates of the total displacement (Fig. F.2). In particular, simulations with either 


< 0.01 or < 0.01 tended to be associated with higher d values (Fig. F.2). The 
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Frequency Frequency Frequency 


estimates of d ranged widely in values, but appeared to by generally higher. 



Figure F.l: Histograms of the parameter estimates for the set of 500 simulations 
of the polar bear sea ice model. The vertical lines represent the parameter values 
used in the simulations. The black lines represent the value used to simulate the 
data ("see Table 


F.l). 
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RMSE with 6 




RMSE with 0 



Figure F.2: Histograms of the RMSE of the states for the polar bear sea ice 
model and of the total displacement, d. The left column represents results when 
the model was htted using the estimated parameter values. The right column 
represents the results when the model was htted using the simulation values. 
The purple columns represent the simulations for which either < 0.01 or 
ae,v < 0 . 01 . 
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G Results when the measurement error is fixed 


As explained in section 2.2 of the main text, we investigated whether fixing the measnrement 
error resolved the parameter estimation problem. To do so, we fitted onr simple likelihood 
(eqn. to the same simnlations as in section 2.1, bnt we fixed the standard deviation of the 
measurement equation to the value used to simulate the data, = 0.1. We only estimated 
the remaining parameters, = (p, <7,,). As for the main analysis, we investigated the 
parameter estimates, RMSE of the states, and likelihood profiles. In addition, we explored 
the likelihood surfaces. 


As mentionned in section 3^ of the main text, fixing the standard deviation of the mea¬ 
surement error to the simulated value, = 0.1, helped reduce the estimation problems 


(compare Fig. G.l to Fig. A.l). In particular, the state RMSE when the parameters were 


estimated were much closer to those when the parameters were fixed to the simulated values 


(e.g., compare Fig. G.ID to Fig. A.ID). In this case, only 5.0% of the simulations had a 


RMSEg value that was 50% larger than their RMSE^. In addition, likelihood profiles were 


more unimodal than when all parameters were estimated (e.g., compare Fig. G.2J to Fig. 
ii). and the state estimates were no longer simply echoing the observations (e.g., compare 
Fig. 1^ to Fig. R)' However, using measurement error information did not completely 
resolve the estimation problems. Some parameter estimates continued to be on the boundary 


of parameter space and far from their simulated values (e.g.. Fig. G.IE). In addition, some 


likelihood profiles remained fiat and some GIs spanned the entire parameter space (e.g.. Fig. 


G.2B). 
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Figure G.l: Histograms of the parameter estimates and of the RMSE of the 
states when the values of the measurement error is hxed to its true value. Each 
row represents the results of 200 simulations for a set of parameter values. For 
the hrst three columns, the vertical lines represent the parameter values used 
in the simulations, with black lines used for the values that remain constant, 
(Te = 0.1 and p = 0.7, and the red lines for values that change between set, 
(Tr, = (0.01,0.02,0.05,0.1,0.2,0.5,1). In the last column, the grey histograms 
represent the RMSE of the model htted using the estimated parameter values, 
while the blue histograms represent the RMSE when the model was fitted using 
the simulation values. 
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Figure G.2: Log likelihood surface and profile for the problematic simulations 
when the standard deviation of the measurement equation is hxed to the simulated 
value, (Je = 0.1. The hrst column represents the log likelihood surface for the 
two estimated parameters, = (p, CTn). The grey dot is the maximum likelihood 
estimate and the grey cross is the simulated value. The second and third columns 
represent the log likelihood prohle for p and The curve represents the log 
likelihood when the focal parameter is hxed (the other parameter are optimise to 
maximise the log likelihood). The dash lines are the true parameter values (i.e., 
value used for the simulation), the full lines are the maximum likelihood estimates 
and the grey bands represent the 95% Cl. The last column shows the time-series. 
The black lines represent the observations, yt, the red lines the simulated true 
states, Xt, and the grey dashed lines the estimated states, Xt- 
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H Diagnostic tools 


Identifying whether the model has parameter estimability problems is an important step 
towards avoiding biased inference. In this Appendix, we present a few diagnostic tools for 
both the more traditional likelihood-based methods and for Bayesian methods. 


One of the best way to check whether a model is capable of estimating the parameters and 
states is throngh a simnlation stndy snch as the one presented in this manuscript. While 
an extensive simulation study that investigates a variety of parameter values is necessary 
to assess the overall capacity of the model, in many cases it may be sufficient to focus on 
parameter values similar to those estimated from the real data. An example of such focussed 
simulation study is presented in Appendix One important aspect of these simulation 
studies is that they require to run a large sample of simulations (we used a sample of 200- 
500 simulations per parameter set). Using a single simulation can be extremely misleading. 


For instance, only 29.6% of our main simulations were problematic (see section 3.1 from the 
main text). However, repeatedly fitting a model to multiple time-series may only be feasible 
with computationally-efficient method. Computing a simulation study with rjags is much 
less practical than with TMB. 


As shown repeatedly in the manuscript, one way to identify the potential for estimation prob¬ 
lems with likelihood-based methods is by investigating the profile likelihood. Flat, jagged, 
or bimodal profile likelihoods indicate the potential for parameter-estimation problems (e.g.. 
Fig. [^-C). In contrast, smooth unimodal profile likelihoods, such as those where we know 
the true states (e.g.. Fig. |D.1[ 4-C), indicate that there is no obvious estimation problems. 
A more comprehensive investigation can be done through the visualization of a likelihood 


surface (Polansky et ah, 2009). For example, by looking at the parameters two-by-two. As 
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an example, we used one of the problematic simulations (i.e., RMSE^ > 1.5 x RMSE^), 
where arj = 0.05. We computed the likelihood surface for the measurement error and pro¬ 
cess stochasticity. To demonstrate the difference between a problematic and a well-behaved 
model, we compared the case when the states were estimated (as in the main text) to the 
case when the states were known (Appendix]^. We can see that when the states are esti¬ 
mated the likelihood surface has a diagonal ridge indicating that it is difficult to separate 
the values of cr^ from those of (Fig. |H.1| A). In contrast, in the case when the states are 
known, we have a well-behaved unimodal likelihood surface (Fig. H.lp). 


When the model is htted with Bayesian methods, chain convergence is often used as a diag¬ 
nostic. The sample paths of MCMC chains for non-identihable parameters may interchange 


their values and lead to numerical or convergence problems (Cressie et ah, 2009). However, 


in our case, very few replicates had convergence problems and many of the converged chains 
lead to biased estimates (see Appendix]^. To further investigate the potential for estimation 
problems, in particular, to verify whether the estimates from the different parameters are 
correlated, one can investigate the posterior distribution of parameters. To show how this 
method has similarity to investigating the likelihood surface, we used the same example as 
above. We compared the posterior distribution of model described in eq. [I]|^ (see Appendix 
[C| for the description of priors) when the states were estimated as opposed to when the 


states were known. As we can see in Fig. H.2 4, when the states are estimated the posterior 
distribution of the measurement error and process stochasticity appears strongly correlated, 
indicating that there is an estimation problem. In contrast, the posterior distribution when 
the states are known does not appear correlated (Fig. H.2p ), indicating that there is no obvi¬ 
ous estimation problem. This is consistent with the results of Appendix which shows that 














when we know the true states, the estimated parameters are close to their true values. 



Figure H.l: Log likelihood surface for the measurement error and process stochas- 
ticity. A) Surface when the states and parameters are estimated. B) Surface when 
the true states are known. 




Figure H.2: Posterior distribution of the measurement error and process stochas- 
ticity. A) Distribution when the states and parameters are estimated. B) Distri¬ 
bution when the true states are known. 
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